Transient behavior of heat transport in a thermal switch 
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We study the time-dependent transport of heat in a nanoscale thermal switch. The switch consists 
of left and right leads that are initially uncoupled. During switch-on the coupling between the leads is 
abruptly turned on. We use the nonequilibrium Green's function formalism and numerically solve the 
constructed Dyson equation to determine the nonperturbative heat current. At the transient regime 
we find that the current initially flows simultaneously into both of the leads and then afterwards 
oscillates between flowing into and out of the leads. At later times the oscillations decay away and 
the current settles into flowing from the hotter to the colder lead. We find the transient behavior 
to be influenced by the extra energy added during switch-on. Such a transient behavior also exists 
even when there is no temperature difference between the leads. The current at the long-time limit 
approaches the steady-state value independently calculated from the Landauer formula. 

PACS numbers: 44.10.+i,63.22.-m,66.70.-f,66.70.Lm 
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The physics of generation, dissipation, and manipula- 
tion of heat in nanoscale systems is an important topic 
that has recently gathered attention. Experiments on 
molecular junctions" found that the heat generated in 
current-carrying metal-molecule junctions is substantial 
and can affect the integrity of the device. Understanding 
how to efficiently dissipate extraneous heat in nanoscale 
systems is thus imperative in the construction of devices. 
In addition, the movement of heat may be harnessed for 
information processing^. Recent experiments have shown 
the viability of thermal transistors 3 , thermal rectifiers us- 
ing inhomogeneous carbon and boron nitride nanotubes 4 , 
and conductance-tunable thermal links consisting of mul- 
tiwalled carbon nanotubes^. 

Most of the above-mentioned work, however, focus on 
the examination of steady-state phenomena. In contrast, 
any physical device must function in a time-dependent 
environment. Although some work has been done in the 
study of time-dependent electronic transport^, the in- 
vestigation of time-dependent behavior in quantum heat 
transport has not yet attracted much attention!. De- 
veloping theoretical tools and computational methods 
for such problems are thus essential to the progress of 
the field. In this paper we study the time-dependent 
heat current in a junction system we call a thermal 
switch. We generalize the nonequilibrium Green's func- 
tion formalism 8 , which is developed for steady-state sit- 
uations, to the time-dependent case with a well-defined 
initial thermal state. To get nonperturbative results, the 
key steps we take are to construct and numerically solve 
a Dyson equation that does not satisfy time-translational 
invariance. 

Fig.[T]shows a one-dimensional chain having a coupling 
that can be switched on and off. The semi-infinite left 
and right leads are linear chains of masses m. Each atom 
interacts with its left and right neighbors through an in- 
terparticle harmonic potential having spring constant k. 
An on-site harmonic potential, with spring constant fco, is 
also experienced by each atom. During time t < the left 
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FIG. 1: (color online) An illustration of a quantum thermal 
switch. The labels of the first 3 sites in each lead are shown. 
The coupling between sites and 1 is switched on at t — 0. 



and right leads are uncoupled and are in thermal equilib- 
rium with temperatures TL and Tr, respectively. At time 
t = the coupling potential, in the form of an interpar- 
ticle harmonic potential with the same spring constant 
k, is switched on, i.e., the potential between the masses 
labeled and 1 in Fig. [1] is suddenly switched on. We 
then want to know how the time-dependent heat current 
behaves in such a setup. In experiments on molecular 
junctions^ a scanning tunneling microscope tip is used to 
stretch a molecule until a bond in the molecule breaks. In 
the thermal switch we can think of the switch-on as the 
inverse process, i.e., a bond is induced through proximity. 
The leads follow the Hamiltonian 



H' 



ij 



(1) 



where the sums are over all sites in the lead, the trans- 
formed coordinates are given by Ui = ^/mxi, Xi is the 
relative displacement of the i-th atom of mass m, and 
K a is the spring constant matrix. The K L and K K ma- 
trices are semi-infinite tridiagonal matrices with 2k + fco 
along the diagonal and — k along the off-diagonal. The 
Hamiltonian for the switched coupling is 



H LR = 



E 



(2) 



The coupling constant matrices V hK and V Rh are zero 
matrices during t < 0. After the switch-on, V LK has 
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one non-zero element V^ R = — k and V RL has the lone 
non-zero element = —k, where the matrix indices 

correspond to the labels of the masses in the leads. Note 
that H hR — H RL . The time-dependent governing Hamil- 
tonian therefore is H(t) = H L + H R + H LK 9 (t), where 
9 (t) is the Heaviside step function. 

The current flowing out of the left lead I^(t) = 
— (dH^/dt), i.e., it is the expectation value of the rate of 
change in iJ L . When the switch is turned on the current 
is given by 



J L (i) = hk Im 



dG^<{hM) 
dt 2 



(3) 



where "Im" stands for the imaginary part. The lesser 
Green's function that appears in the formula is defined 
as G RL ' < (t 1 ,t 2 ) = -j-{u^ (t 2 )uf (h)), where the sub- 
scripts and 1 correspond to the labels of the masses. 
Note that this is a two-time correlation function that does 
not satisfy time-translational invariance, i.e., its time- 
dependence can not be written as a difference t\ — t 2 . 
Similarly, the current flowing out of the right lead, 
In(t) = — (dH n /dt\, when the switch is turned on is in 
the form of Eq. ^ except that the R and L superscripts 
are swapped. 

To determine the full, nonperturbative, current we are 
going to solve the associated Dyson equation. First, we 
define the contour-ordered Green's function^ 



G RL (r 1 ,r 2 ) = --<T c uf (nK (7a)) 



(4) 



where the u R and Uq are Heisenberg operators, T c is the 
contour-ordering operator, and t\ and t 2 are complex 
variables on the contour C . To determine the current at 
time t we employ a Keldysh contour C that goes from 
time to t and then back to 0. Since the left and right 
leads are uncorrelated before the switch is turned on at 
time 0, the contour does not have a complex tail after it 
goes back to time 0. 

Converting to the interaction picture the contour- 
ordered Green's function shown in Eq. (U) becomes 



t/ C ^ LR ^\H (Ti)u L (T2) 



G^{r x ,r 2 ) = ~{T c e 

(5) 

A perturbative calculation can then be done by expand- 
ing the exponential as an infinite series. We can also use 
this expansion in constructing the Dyson equation. In 
the series, the zeroth-order term vanishes because it does 
not contain the coupling potential. Furthermore, all the 
even-ordered terms also vanish because there will be an 
extra m r -m l pair without a connecting coupling poten- 
tial. Thus, only the odd-ordcrcd terms survive. 

We can use diagrams in constructing the Dyson equa- 
tion. Shown in Fig. [2] is the resulting diagram equa- 
tion when we expand the series in Eq. ([3]). A double- 
line diagram represents G RL , a single line represents 
the equilibrium Green's functions for the right lead, 
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FIG. 2: (color online) Diagram representation of the Dyson 
equation. Each line is labelled by the Green's function it 
represents. Each concentric dot represents a coupling vertex. 



g R {Ti,T 2 ) = -j- (T c uf (ri)uf (t 2 )) q , and a dashed line 
represents the equilibrium Green's functions for the left 
lead, <7 L (ti,t 2 ) = (T c u L (n) u L (r 2 )) . The sub- 
script implies that the average is taken with respect 
to equilibrium distributions that are maintained when 
t < before the switch-on. 

Rewriting the diagram equation in Fig. [2] as contour 
integrals, we have 

G RL (n,r 2 ) = / dr a g R (ti, To) V RL g L (r„, r 2 ) 
Jc 

+ dr a dn g R ( Tl ,T a )V 
Jc Jc 



RL 



IC JC 



x g L (T a ,r b )V LR G RL (r b ,T 2 ). (6) 

Applying Langreth's theorem to Eq. © and then 
iterating^, we should obtain an expression for G RL,< . To 
calculate the current in Eq. ([3]) we need the time deriva- 
tive of G RL:< , and so we differentiate it to get 

dG^<(h,t 2 ) _ dGf h '< (h,t 2 ) 



dt 2 



dt 2 



RL.r U t \ 9Gf L,< (t a ,t 2 ) 



dU 



- k [ dt a G^' 1 (ti,t a ) 

Jo 

- kfdt a Gf-<{t x ,t a ) gg^^g) 

Jo ot 2 

+ k 2 [ dt a f dt b G RL ' r (fi,t a ) 
Jo Jo 

x G^<(t a ,t b ) dGRh Z {tb,t2 \ (7) 



dt 2 



where 



G R (h,t 2 ) = -k / dt a { 5 Rr {h - t a ) g L '< (*« - ta) 



+ g^ih-ta) g L -*(t a -t 2 )} (8) 



is the first-order term in the perturbation series in 
Eq. ([5]). The analytic expressions for the equilibrium 
surface Green's functions <? Rr , g R,< , ff L,< , and <? L ' a are 
known in frequency spaced. To determine their time- 
dependence we numerically calculate their corresponding 
Fourier transforms. 
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The other unknowns in Eq. (J7J involve the retarded 
and advanced versions of the full Green's function. We 
can apply Langreth's theorem again to Eq. (j6|) to deter- 
mine expressions for these unknowns. We get 

G^(h,t 2 ) = -k f dt a G^ p (h,t a )G^(t a ,t 2 ) 
Jo 

+ Gf h ' (h,t 2 ) , (9) 
where ft — r, a, and the first-order term is 

Gf L ^ (h,t 2 ) = -k f dta g^ih-ta) g L ^(t a -t 2 ). 
Jo 

(10) 

To solve Eq. © we discretize the time variable into N 
segments and thus transforming the integral into a sum. 
This results in a linear problem of the form Ax = b, 
where the unknown x is determined by performing an 
LU decomposition on A and then using b in the back- 
substitution. The <9G RL,a (t a , t 2 ) /dt 2 term required in 
Eq. ([7]) can also be calculated by first differentiating 
Eq. © and then finding the solution to the resulting 
equation with the time discretized. By determining the 
time derivative of the full Green's function in Eq. ([7]) 
the current that we calculate is a nonperturbative result. 
We follow the same steps to independently calculate the 
current flowing out of the right lead. 

Shown in Fig. |3] are plots of the time-dependent cur- 
rent flowing out of the leads when the left and right lead 
temperatures are Tl = 330 K and Tr, = 270 K, respec- 
tively. The average temperature is thus T = 300 K 
with offsets of ±10%. The currents oscillate at a fre- 
quency comparable to the highest phonon frequencies 
available in the system and then gradually decay to 



their steady-state values. The dots shown at the right 
edges of the plots are steady-state values calculated in- 
dependently from the Landauer formula 7l = —Ik — 
h J-col% huj (/ L ~ /iO^M* where f L and f R are the 
Bose-Einstein distributions of the left and right leads, 
respectively, and 6 (u>) is 1 within the phonon band, 
ko < oj 2 < 4k + ko, and otherwise^. At steady-state, 
heat should flow from the hotter to the colder lead, i.e., 
from the left to the right lead. Thus, the sign of the cur- 
rent flowing out of the left lead should be positive and 
for the right lead negative. However, during the tran- 
sient time, the current can flow in unexpected directions. 
Just after the switch-on, the current actually does not 
flow from the hotter to the colder lead. In Fig. [3] we 
see the current to flow simultaneously into both of the 
leads. There appears to be an energy source in-between 
the leads that supply the current. Recall that there is 
no coupling between the leads before the switch-on. By 
turning on the switch we actually add energy, in the form 
of the switched coupling potential, to the system. This 
added coupling energy supplies the current that flows 
into both leads. In Fig. [3] we also compare the results 
from first-order perturbation to the results from nonper- 
turbative calculations. Unlike the perturbative results, 
nonperturbative results approach the steady-state values 
at later times. Since the switched coupling has the same 
strength as the interparticle potential we indeed expect 
corrections to perturbative calculations to be significant. 

Shown in Fig. [4] are plots of the sum of the currents, 
Is = Jl + Jr., as functions of the time and the average 
temperature of the leads. Is fluctuates around zero with 
a fluctuation amplitude that decays with time. Fig. @|b) 
shows Is to vary strongly with temperature at the tran- 
sient regime. As time goes on Is slowly loses its de- 
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FIG. 3: The time-dependent current flowing out of the (a) 
left lead and (b) right lead. The (blue online) data points are 
results from solving the Dyson equation while the (red online) 
line is the result from the first-order perturbation calculation. 
The average temperature between the leads is T — 300 K. 
The interparticle spring constant is k = 0.625 eV/(A 2 u) while 
the on-site spring constant is ko = 0.0625 eV/(A 2 u). 



FIG. 4: (a) Plots of Is = II + In as functions of time when 
T = 10 K (orange online), T = 300 K (red online), and 
T = 1000 K (blue online). The temperature offsets are ±10%. 
(b) Plots of |7s| as functions of T at time t = 0.8 [t] (brown 
online), t = 4.1 [t] (yellow online), t = 6.0 [t] (violet online), 
and t = 23.0 [t] (green online), where [t] = 10 -14 s. (c) An 
enlarged view of (a) for time t — 20 [t] to t — 40 [t]. 



4 



60 




-120 M 

-140 I 1 1 ' 

10 20 30 40 

time [10" 14 s] 

FIG. 5: The current flowing through the junction when the 
left and right leads have the same temperature. T = 10 K 
(red online) and T = 300 K (blue online). 

pendence on temperature. Fig. Hfc) shows that at later 
times fluctuations in Is still appear but are significantly 
smaller than those at the transient regime. We do expect 
that at the steady state, of which the long-time limit of 
our data approaches, all of the current flowing out of the 
left lead should flow into the right lead and thus result- 
ing in 7g = 0, regardless of the value of the temperature. 
However, since energy is added during the switch-on, this 
extra energy influences the transient behavior of the sys- 
tem until it eventually dissipates out to the heat baths. 
Fitting the envelope to an exponential function we find a 
characteristic decay time of about 10 x 10 -14 s. At any 
particular time t we can determine the energy the system 
absorbs or emits from (i/ LR (£)) = J Q (Jl, + Ir) dt. 

Suppose we set the temperatures of the leads to be 
the same. Shown in Fig. [S] is the time-dependent cur- 
rent, in either lead since the leads are indistinguish- 
able, for such a situation. At the steady state there 
should be no current flowing within the system. How- 



ever, since energy is added to the system during the 
switch-on, at the transient regime we see a fluctuating 
current with temperature-dependent amplitude flowing 
within the system. 

The method can also be generalized in a straightfor- 
ward manner to deal with a time- varying coupling that 
is on during t > 0. For a mildly increasing coupling such 
as, for example, k(t) = fctanh(/i), the transient current 
also initially flows simultaneously into both of the leads. 
This transient behavior persists because switching on the 
coupling introduces energy, that has to dissipate into the 
baths, into the system. 

To summarize, we have shown an exact nonperturba- 
tive method to calculate the time-dependent heat current 
in a thermal switch using nonequilibrium Green's func- 
tions. The Dyson equation is constructed using Keldysh 
contours and the real-time Green's functions needed to 
calculate the current are determined by applying Lan- 
greth's theorem to the Dyson equation and numerically 
solving the equation with the discretized time variable. 
We set the strength of the switched coupling to be the 
same as the interparticle spring constant. Nonperturba- 
tive results are thus significantly different from pertur- 
bative ones. We find the transient current just after the 
switch-on to be influenced by the extra switched coupling 
energy. In particular, the initial reaction after switch-on 
is for the current to flow simultaneously into both the 
left and right leads. The current then oscillates with am- 
plitude that decays with time. In the long-time limit the 
current approaches the expected steady-state values cal- 
culated independently from the Landauer formula. We 
also note that the theory presented here is not restricted 
only to one-dimensional chains but is applicable to any 
junction system where a thermal switch-on occurs. 
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